pacman::p_load(tidyverse, patchwork, ggcorrplot)
path1      <-
    '~/Git/henriquelaureano.github.io/THESIS/code/coefs/coefs1to100/'
path2      <-
    '~/Git/henriquelaureano.github.io/THESIS/code/coefs/coefs101to200/'
fixednames <- c('beta1', 'beta2', 'gama1', 'gama2', 'w1', 'w2')
varnames   <- c('logs2_1', 'logs2_2', 'logs2_3', 'logs2_4')
cornames   <- c(
    'rhoZ12', 'rhoZ13', 'rhoZ14', 'rhoZ23', 'rhoZ24', 'rhoZ34'
)
metnames   <- c('conv', 'mll')

vdata <- function(data, path)
{
    out <-
        tibble::as_tibble(read.table(paste0(path, data)))%>%
        dplyr::select(-V1)
    type <- substring(data, 6, 7)    
    switch(type,
           v1={ out <-
                    out%>%'colnames<-'(c(
                              fixednames,
                              varnames[1:2], cornames[1], metnames)) },
           v2={ out <-
                    out%>%'colnames<-'(c(
                              fixednames,
                              varnames[3:4], cornames[6], metnames)) },
           v3={ out <-
                    out%>%'colnames<-'(c(
                              fixednames,
                              varnames, cornames[c(1, 6)], metnames)) },
           v4={ out <-
                    out%>%'colnames<-'(c(
                              fixednames,
                              varnames, cornames, metnames)) },
           stop('Invalid type'))
    return(out)
}

v1data <- v2data <- v3data <- v4data <- vector('list', 18)
label  <- numeric(18)
l      <- 1
for (i in 1:3)
{
    for (j in 1:2)
    {
        for (k in 1:3)
        {
            label[l]    <- paste0('cs', i, '_', j, '_', k)
            v1data[[l]] <-
                dplyr::bind_rows(
                           vdata(paste0('coefsv1_', label[l], '.txt'),
                                 path1),
                           vdata(paste0('coefsv1_', label[l], '.txt'),
                                 path2)
                       )
            v2data[[l]] <-
                dplyr::bind_rows(
                           vdata(paste0('coefsv2_', label[l], '.txt'),
                                 path1),
                           vdata(paste0('coefsv2_', label[l], '.txt'),
                                 path2)
                       )
            v3data[[l]] <- vdata(paste0('coefsv3_', label[l], '.txt'),
                                 path1)
            v4data[[l]] <- vdata(paste0('coefsv4_', label[l], '.txt'),
                                 path1)
            l <- l+1
        }
    }
}
## str(v1data);length(v1data)
## str(v2data);length(v2data)
## str(v3data);length(v3data)
## str(v4data);length(v4data)

cif1    <- setNames(c( 3,  2.6, 2.5, 4, 5, 10), fixednames)
cif2    <- setNames(c(-2, -1.5, 1, 1.5, 3, 4), fixednames)
logvars <- setNames(log(c(1, 0.6, 0.7, 0.9)), varnames)
rhoZs   <- setNames(atanh(c(0.1, -0.5, 0.3, 0.3, -0.4, 0.2)), cornames)

vtrue <- function(data.cif1, data.cif2)
{
    out <- rep(c(replicate(data.cif1, n=3, simplify=FALSE),
                 replicate(data.cif2, n=3, simplify=FALSE)), 3)
    return(out)
}
v1cif1 <- c(cif1, logvars[1:2], rhoZs[1])
v1cif2 <- c(cif2, logvars[1:2], rhoZs[1])
v1true <- vtrue(v1cif1, v1cif2)
v2cif1 <- c(cif1, logvars[3:4], rhoZs[6])
v2cif2 <- c(cif2, logvars[3:4], rhoZs[6])
v2true <- vtrue(v2cif1, v2cif2)
v3cif1 <- c(cif1, logvars, rhoZs[c(1, 6)])
v3cif2 <- c(cif2, logvars, rhoZs[c(1, 6)])
v3true <- vtrue(v3cif1, v3cif2)
v4cif1 <- c(cif1, logvars, rhoZs)
v4cif2 <- c(cif2, logvars, rhoZs)
v4true <- vtrue(v4cif1, v4cif2)
cerror <- function(data, true, data_label) {
    coefs <- data%>%dplyr::select(-c(conv, mll))
    error <- coefs
    for (i in seq(ncol(coefs))) { error[i] <- true[i]-coefs[i] }
    out <- error%>%
        dplyr::summarize_all(mean)%>%
        dplyr::add_row(
                   error%>%
                   dplyr::summarize_all(quantile, c(0.025, 0.975)))%>%
        dplyr::add_row(
                   error%>%
                   dplyr::summarize_all(sd))%>%
        dplyr::mutate(label=c('mean', 'q025', 'q975', 'sd'),
                      conv=1-mean(data$conv),
                      nr=nrow(data), 
                      data=data_label)
    return(out)
}

label <- sub('cs1_', 'cs02-',  label)
label <- sub('cs2_', 'cs05-',  label)
label <- sub('cs3_', 'cs10-', label)

label <- sub('1_', 'high-', label)
label <- sub('2_', '-low-',  label)

label <- sub('-1$', '-05k',  label)
label <- sub('-2$', '-30k', label)
label <- sub('-3$', '-60k', label)

v1error <- v2error <- v3error <- v4error <- NULL
for (i in seq(18))
{
    v1error <-
        v1error%>%
        dplyr::bind_rows(cerror(v1data[[i]], v1true[[i]], label[i]))
    v2error <-
        v2error%>%
        dplyr::bind_rows(cerror(v2data[[i]], v2true[[i]], label[i]))
    v3error <-
        v3error%>%
        dplyr::bind_rows(cerror(v3data[[i]], v3true[[i]], label[i]))
    v4error <-
        v4error%>%
        dplyr::bind_rows(cerror(v4data[[i]], v4true[[i]], label[i]))
}

## v1error
## v2error
## v3error
## v4error

inside <- function(par, errordata)
{
    out <-
        errordata%>%
        dplyr::select(par, label, conv, nr, data)%>%
        tidyr::pivot_wider(names_from=label, values_from=par)
    out$data <-
        forcats::fct_relevel(forcats::as_factor(out$data), rev)
    return(out)
}
biasformat <- function(par, errordatas)
{
    out <- dplyr::bind_rows(inside(par, errordatas[[1]]),
                            inside(par, errordatas[[2]]),
                            inside(par, errordatas[[3]]),
                            inside(par, errordatas[[4]]),
                            .id='v')
    out$v <- forcats::fct_recode(out$v,
                                 'RISK MODEL'      ='1',
                                 'TIME MODEL'      ='2',
                                 'BLOCK-DIAG MODEL'='3',
                                 'COMPLETE MODEL'  ='4')
    return(out)
}

bias2plot <- function(df, title)
{
    ggplot(df, aes(ymin=q025, ymax=q975, x=data, color=conv, size=nr))+
        geom_linerange(position=position_dodge(width=0.2))+
        geom_point(mapping=aes(x=data, y=mean), size=3,
                   shape=21, color='black', fill='white', stroke=2)+
        facet_grid(~v, scales='free_x')+
        geom_hline(yintercept=0, linetype='dashed')+
        labs(x=NULL, y='Bias',
             title=title,
             subtitle='with 2.5% and 97.5% quantiles',
             color='Convergance\n(%)',
             size='Number of\nmodels')+
        coord_flip()+
        theme_minimal()+
        scale_color_viridis_c()+
        ## scale_color_distiller(palette='Greys', direction=1)+
        guides(color=guide_colorbar(barwidth=20))+
        theme(
            strip.background=element_rect(colour='black', fill='white'),
            legend.position='bottom',
            plot.title=element_text(size=15), 
            plot.subtitle=element_text(size=14, face='italic'),
            axis.text.x=element_text(size=12), 
            axis.text.y=element_text(size=13),
            axis.title.x=element_text(size=13,
                                      margin=unit(c(t=3, r=0, b=0, l=0),
                                                  'mm')),
            legend.justification=c(1, 0),
            legend.title=element_text(size=13), 
            legend.text=element_text(size=12),
            strip.text.x=element_text(size=13, face='bold')
        )
}

bias <- function(par, title)
{
    out <- biasformat(
        par,
        list(v1error, v2error, v3error, v4error)
    )
    bias2plot(out, title=title)
}

bias('beta1', expression(bold('Parameter:'~beta[1])))

bias('beta2', expression(bold('Parameter:'~beta[2])))

bias('gama1', expression(bold('Parameter:'~gamma[1])))

bias('gama2', expression(bold('Parameter:'~gamma[2])))

bias('w1', expression(bold('Parameter:'~w[1])))

bias('w2', expression(bold('Parameter:'~w[2])))

par <- 'logs2_1'
out <- dplyr::bind_rows(inside(par, v1error),
                        inside(par, v3error),
                        inside(par, v4error),
                        .id='v')
out$v <- forcats::fct_recode(out$v,
                             'RISK MODEL'      ='1',
                             'BLOCK-DIAG MODEL'='2',
                             'COMPLETE MODEL'  ='3')
bias2plot(out, expression(bold('Parameter:'~log(sigma[1]^2))))

par <- 'logs2_2'
out <- dplyr::bind_rows(inside(par, v1error),
                        inside(par, v3error),
                        inside(par, v4error),
                        .id='v')
out$v <- forcats::fct_recode(out$v,
                             'RISK MODEL'      ='1',
                             'BLOCK-DIAG MODEL'='2',
                             'COMPLETE MODEL'  ='3')
bias2plot(out, expression(bold('Parameter:'~log(sigma[2]^2))))

par <- 'logs2_3'
out <- dplyr::bind_rows(inside(par, v2error),
                        inside(par, v3error),
                        inside(par, v4error),
                        .id='v')
out$v <- forcats::fct_recode(out$v,
                             'TIME MODEL'      ='1',
                             'BLOCK-DIAG MODEL'='2',
                             'COMPLETE MODEL'  ='3')
bias2plot(out, expression(bold('Parameter:'~log(sigma[3]^2))))

par <- 'logs2_4'
out <- dplyr::bind_rows(inside(par, v2error),
                        inside(par, v3error),
                        inside(par, v4error),
                        .id='v')
out$v <- forcats::fct_recode(out$v,
                             'TIME MODEL'      ='1',
                             'BLOCK-DIAG MODEL'='2',
                             'COMPLETE MODEL'  ='3')
bias2plot(out, expression(bold('Parameter:'~log(sigma[4]^2))))

par <- 'rhoZ12'
out <- dplyr::bind_rows(inside(par, v1error),
                        inside(par, v3error),
                        inside(par, v4error),
                        .id='v')
out$v <- forcats::fct_recode(out$v,
                             'RISK MODEL'      ='1',
                             'BLOCK-DIAG MODEL'='2',
                             'COMPLETE MODEL'  ='3')
bias2plot(out, expression(bold('Parameter:'~z(rho[12]))))

par <- 'rhoZ34'
out <- dplyr::bind_rows(inside(par, v2error),
                        inside(par, v3error),
                        inside(par, v4error),
                        .id='v')
out$v <- forcats::fct_recode(out$v,
                             'TIME MODEL'      ='1',
                             'BLOCK-DIAG MODEL'='2',
                             'COMPLETE MODEL'  ='3')
bias2plot(out, expression(bold('Parameter:'~z(rho[34]))))

bias2plot_rho <- function(df, title)
{
    ggplot(df, aes(ymin=q025, ymax=q975, x=data, color=conv, size=nr))+
        geom_linerange(position=position_dodge(width=0.2))+
        geom_point(mapping=aes(x=data, y=mean), size=3,
                   shape=21, color='black', fill='white', stroke=2)+
        facet_grid(~v, scales='free_x', labeller=label_parsed)+
        geom_hline(yintercept=0, linetype='dashed')+
        labs(x=NULL, y='Bias',
             title=title,
             subtitle='with 2.5% and 97.5% quantiles',
             color='Convergance\n(%)',
             size='Number of\nmodels')+
        coord_flip()+
        theme_minimal()+
        scale_color_viridis_c()+
        guides(color=guide_colorbar(barwidth=12.5))+
        theme(
            strip.background=element_rect(colour='black', fill='white'),
            legend.position='bottom',
            plot.title=element_text(size=15), 
            plot.subtitle=element_text(size=14, face='italic'),
            axis.text.x=element_text(size=12), 
            axis.text.y=element_text(size=13),
            axis.title.x=element_text(size=13,
                                      margin=unit(c(t=3, r=0, b=0, l=0),
                                                  'mm')),
            legend.justification=c(1, 0),
            legend.title=element_text(size=13), 
            legend.text=element_text(size=12),
            strip.text.x=element_text(size=13)
        )
}

out <- dplyr::bind_rows(inside('rhoZ13', v4error),
                        inside('rhoZ24', v4error),
                        inside('rhoZ14', v4error),
                        inside('rhoZ23', v4error),
                        .id='v')
out$v <- forcats::fct_recode(out$v,
                             'bold(z(rho[13]))'='1',
                             'bold(z(rho[24]))'='2',
                             'bold(z(rho[14]))'='3',
                             'bold(z(rho[23]))'='4')
bias2plot_rho(out,
              expression(bold("Complete model's cross-correlations")))

bias2plot2 <- function(df, title)
{
    ggplot(df, aes(ymin=mean-1.96*sd, ymax=mean+1.96*sd, x=data,
                   color=conv, size=nr))+
        geom_linerange(position=position_dodge(width=0.2))+
        geom_point(mapping=aes(x=data, y=mean), size=3,
                   shape=21, color='black', fill='white', stroke=2)+
        facet_grid(~v, scales='free_x')+
        geom_hline(yintercept=0, linetype='dashed')+
        labs(x=NULL, y='Bias',
             title=title,
             subtitle='with \u00b1 1.96 standard deviations',
             color='Convergance\n(%)',
             size='Number of\nmodels')+
        coord_flip()+
        theme_minimal()+
        scale_color_viridis_c()+
        ## scale_color_distiller(palette='Greys', direction=1)+
        guides(color=guide_colorbar(barwidth=20))+
        theme(
            strip.background=element_rect(colour='black', fill='white'),
            legend.position='bottom',
            plot.title=element_text(size=15), 
            plot.subtitle=element_text(size=14, face='italic'),
            axis.text.x=element_text(size=12), 
            axis.text.y=element_text(size=13),
            axis.title.x=element_text(size=13,
                                      margin=unit(c(t=3, r=0, b=0, l=0),
                                                  'mm')),
            legend.justification=c(1, 0),
            legend.title=element_text(size=13), 
            legend.text=element_text(size=12),
            strip.text.x=element_text(size=13, face='bold')
        )
}
bias2 <- function(par, title)
{
    out <- biasformat(
        par,
        list(v1error, v2error, v3error, v4error)
    )
    bias2plot2(out, title=title)
}
bias2('beta1', expression(bold('Parameter:'~beta[1])))

bias2('beta2', expression(bold('Parameter:'~beta[2])))

bias2('gama1', expression(bold('Parameter:'~gamma[1])))

bias2('gama2', expression(bold('Parameter:'~gamma[2])))

bias2('w1', expression(bold('Parameter:'~w[1])))

bias2('w2', expression(bold('Parameter:'~w[2])))

par <- 'logs2_1'
out <- dplyr::bind_rows(inside(par, v1error),
                        inside(par, v3error),
                        inside(par, v4error),
                        .id='v')
out$v <- forcats::fct_recode(out$v,
                             'RISK MODEL'      ='1',
                             'BLOCK-DIAG MODEL'='2',
                             'COMPLETE MODEL'  ='3')
bias2plot2(out, expression(bold('Parameter:'~log(sigma[1]^2))))

par <- 'logs2_2'
out <- dplyr::bind_rows(inside(par, v1error),
                        inside(par, v3error),
                        inside(par, v4error),
                        .id='v')
out$v <- forcats::fct_recode(out$v,
                             'RISK MODEL'      ='1',
                             'BLOCK-DIAG MODEL'='2',
                             'COMPLETE MODEL'  ='3')
bias2plot2(out, expression(bold('Parameter:'~log(sigma[2]^2))))

par <- 'logs2_3'
out <- dplyr::bind_rows(inside(par, v2error),
                        inside(par, v3error),
                        inside(par, v4error),
                        .id='v')
out$v <- forcats::fct_recode(out$v,
                             'TIME MODEL'      ='1',
                             'BLOCK-DIAG MODEL'='2',
                             'COMPLETE MODEL'  ='3')
bias2plot2(out, expression(bold('Parameter:'~log(sigma[3]^2))))

par <- 'logs2_4'
out <- dplyr::bind_rows(inside(par, v2error),
                        inside(par, v3error),
                        inside(par, v4error),
                        .id='v')
out$v <- forcats::fct_recode(out$v,
                             'TIME MODEL'      ='1',
                             'BLOCK-DIAG MODEL'='2',
                             'COMPLETE MODEL'  ='3')
bias2plot2(out, expression(bold('Parameter:'~log(sigma[4]^2))))

par <- 'rhoZ12'
out <- dplyr::bind_rows(inside(par, v1error),
                        inside(par, v3error),
                        inside(par, v4error),
                        .id='v')
out$v <- forcats::fct_recode(out$v,
                             'RISK MODEL'      ='1',
                             'BLOCK-DIAG MODEL'='2',
                             'COMPLETE MODEL'  ='3')
bias2plot2(out, expression(bold('Parameter:'~z(rho[12]))))

par <- 'rhoZ34'
out <- dplyr::bind_rows(inside(par, v2error),
                        inside(par, v3error),
                        inside(par, v4error),
                        .id='v')
out$v <- forcats::fct_recode(out$v,
                             'TIME MODEL'      ='1',
                             'BLOCK-DIAG MODEL'='2',
                             'COMPLETE MODEL'  ='3')
bias2plot2(out, expression(bold('Parameter:'~z(rho[34]))))

bias2plot2_rho <- function(df, title)
{
    ggplot(df, aes(ymin=mean-1.96*sd, ymax=mean+1.96*sd, x=data,
                   color=conv, size=nr))+
        geom_linerange(position=position_dodge(width=0.2))+
        geom_point(mapping=aes(x=data, y=mean), size=3,
                   shape=21, color='black', fill='white', stroke=2)+
        facet_grid(~v, scales='free_x', labeller=label_parsed)+
        geom_hline(yintercept=0, linetype='dashed')+
        labs(x=NULL, y='Bias',
             title=title,
             subtitle='with \u00b1 1.96 standard deviations',
             color='Convergance\n(%)',
             size='Number of\nmodels')+
        coord_flip()+
        theme_minimal()+
        scale_color_viridis_c()+
        guides(color=guide_colorbar(barwidth=12.5))+
        theme(
            strip.background=element_rect(colour='black', fill='white'),
            legend.position='bottom',
            plot.title=element_text(size=15), 
            plot.subtitle=element_text(size=14, face='italic'),
            axis.text.x=element_text(size=12), 
            axis.text.y=element_text(size=13),
            axis.title.x=element_text(size=13,
                                      margin=unit(c(t=3, r=0, b=0, l=0),
                                                  'mm')),
            legend.justification=c(1, 0),
            legend.title=element_text(size=13), 
            legend.text=element_text(size=12),
            strip.text.x=element_text(size=13)
        )
}

out <- dplyr::bind_rows(inside('rhoZ13', v4error),
                        inside('rhoZ24', v4error),
                        inside('rhoZ14', v4error),
                        inside('rhoZ23', v4error),
                        .id='v')
out$v <- forcats::fct_recode(out$v,
                             'bold(z(rho[13]))'='1',
                             'bold(z(rho[24]))'='2',
                             'bold(z(rho[14]))'='3',
                             'bold(z(rho[23]))'='4')
bias2plot2_rho(out,
               expression(bold("Complete model's cross-correlations")))

COR

coefs22_cor <- round(cor(coefs22), 1)
coefs36_cor <- round(cor(coefs36), 1)
coefs38_cor <- round(cor(coefs38), 1)
coefs40_cor <- round(cor(coefs40), 1)

cor2plot <- function(cor.data, title, zrhos) {
    ggcorrplot(cor.data, 
               type='lower', 
               lab=TRUE, 
               ## lab_size=5, 
               ## outline.color='white', 
               colors=RColorBrewer::brewer.pal(n=3, name='Greys'),
               title=title,
               p.mat=cor_pmat(cor.data), 
               insig='blank', 
               ggtheme=theme(
                   plot.title=element_text(size=16, face='bold'), 
                   axis.text.x=element_text(size=15), 
                   axis.text.y=element_text(size=15), 
                   axis.title.x=element_blank(),
                   axis.title.y=element_blank(),
                   panel.background=element_blank(),
                   axis.ticks=element_blank(),
                   legend.justification=c(1, 0),
                   legend.position=c(0.4, 0.85),
                   legend.direction='horizontal',
                   legend.title=element_blank(),
                   legend.text=element_text(size=11)))+
        scale_x_discrete(labels=c(expression(beta[2]),
                                  expression(gamma[1]),
                                  expression(gamma[2]),
                                  expression(w[1]),
                                  expression(w[2]),
                                  expression(log(sigma[1]^2)),
                                  expression(log(sigma[2]^2)),
                                  expression(log(sigma[3]^2)),
                                  expression(log(sigma[4]^2)),
                                  zrhos))+
        scale_y_discrete(labels=c(expression(beta[1]), 
                                  expression(beta[2]),
                                  expression(gamma[1]),
                                  expression(gamma[2]),
                                  expression(w[1]),
                                  expression(w[2]),
                                  expression(log(sigma[1]^2)),
                                  expression(log(sigma[2]^2)),
                                  expression(log(sigma[3]^2)),
                                  expression(log(sigma[4]^2)),
                                  head(zrhos, -1)))+
        coord_cartesian()
}

p1 <- cor2plot(coefs22_cor,
               title='model22 parameters correlation',
               zrhos=c(expression(z(rho[12])),
                       expression(z(rho[34])),
                       expression(z(rho[13])),
                       expression(z(rho[24]))))
p2 <- cor2plot(coefs36_cor,
               title='model36 parameters correlation',
               zrhos=c(expression(z(rho[12])),
                       expression(z(rho[34])),
                       expression(z(rho[13])),
                       expression(z(rho[24])),
                       expression(z(rho[14])),
                       expression(z(rho[23]))))
p3 <- cor2plot(coefs38_cor,
               title='model38 parameters correlation',
               zrhos=c(expression(z(rho[12])),
                       expression(z(rho[34])),
                       expression(z(rho[14])),
                       expression(z(rho[23]))))
p4 <- cor2plot(coefs40_cor,
               title='model40 parameters correlation',
               zrhos=c(expression(z(rho[13])),
                       expression(z(rho[24])),
                       expression(z(rho[14])),
                       expression(z(rho[23]))))
(p1|p2)/(p3|p4)

BIAS

true22 <- c(-2, -1.5, 1.2, 1, 3, 5,
            log(0.4), log(0.4), log(0.25), log(0.25),
            atanh(0.15/sqrt(0.4*0.4)),
            atanh(0.1/sqrt(0.25*0.25)),
            atanh(0.05/sqrt(0.4*0.25)),
            atanh(0.05/sqrt(0.4*0.25)))
true36 <- c(-2, -1.5, 1.2, 1, 3, 5,
            log(0.4), log(0.4), log(0.25), log(0.25),
            atanh(0.15/sqrt(0.4*0.4)),
            atanh(0.1/sqrt(0.25*0.25)),
            atanh(0.05/sqrt(0.4*0.25)),
            atanh(0.05/sqrt(0.4*0.25)), 
            atanh(0.2/sqrt(0.4*0.25)),
            atanh(0.2/sqrt(0.4*0.25)))
true38 <- c(-2, -1.5, 1.2, 1, 3, 5,
            log(0.4), log(0.4), log(0.25), log(0.25),
            atanh(0.15/sqrt(0.4*0.4)),
            atanh(0.1/sqrt(0.25*0.25)),
            atanh(0.1/sqrt(0.4*0.25)),
            atanh(0.1/sqrt(0.4*0.25)))
true40 <- c(-2, -1.5, 1.2, 1, 3, 5,
            log(0.4), log(0.4), log(0.25), log(0.25),
            atanh(0.05/sqrt(0.4*0.25)),
            atanh(0.05/sqrt(0.4*0.25)),
            atanh(0.2/sqrt(0.4*0.25)),
            atanh(0.2/sqrt(0.4*0.25)))

cerror <- function(coefs, true) {
    error <- coefs
    for (i in seq(ncol(coefs))) { error[i] <- true[i]-coefs[i] }
    out <- error%>%
        summarize_all(mean)%>%
        add_row(error%>%summarize_all(quantile, c(0.025, 0.975)))%>%
        mutate(label=c('mean', 'q025', 'q975'))%>%
        pivot_longer(!label, names_to='par', values_to='value')%>%
        pivot_wider(names_from=label, values_from=value)
    out$par <- as_factor(out$par)
    return(out)
}
error22 <- cerror(coefs22, true22)
error36 <- cerror(coefs36, true36)
error38 <- cerror(coefs38, true38)
error40 <- cerror(coefs40, true40)

bias2plot <- function(data, title, zrhos){
    ggplot(data, aes(ymin=q025, ymax=q975, x=par))+
        geom_linerange(
            color=RColorBrewer::brewer.pal(n=3, name='Greys')[2],
            position=position_dodge(width=0.2),
            size=3)+
        ylim(c(-1, 1))+
        geom_point(mapping=aes(x=par, y=mean), size=2)+
        geom_hline(yintercept=0, linetype='dashed')+
        labs(x=NULL, y='Bias',
             title=title, subtitle='with 2.5% and 97.5% quantiles')+
        scale_x_discrete(labels=c(expression(beta[1]),
                                  expression(beta[2]),
                                  expression(gamma[1]),
                                  expression(gamma[2]),
                                  expression(w[1]),
                                  expression(w[2]),
                                  expression(log(sigma[1]^2)),
                                  expression(log(sigma[2]^2)),
                                  expression(log(sigma[3]^2)),
                                  expression(log(sigma[4]^2)),
                                  zrhos))+
        coord_flip()+
        theme_minimal()+
        theme(plot.title=element_text(size=12, face='bold'),
              plot.subtitle=element_text(size=11, face='italic'),
              axis.text.x=element_text(size=11),
              axis.text.y=element_text(size=11),
              axis.title.x=element_text(
                  size=12, 
                  margin=unit(c(t=3, r=0, b=0, l=0), 'mm')))
}
p1 <- bias2plot(error22,
                title='model22 parameters bias',
                zrhos=c(expression(z(rho[12])),
                        expression(z(rho[34])),
                        expression(z(rho[13])),
                        expression(z(rho[24]))))
p2 <- bias2plot(error36,
                title='model36 parameters bias',
                zrhos=c(expression(z(rho[12])),
                        expression(z(rho[34])),
                        expression(z(rho[13])),
                        expression(z(rho[24])),
                        expression(z(rho[14])),
                        expression(z(rho[23]))))
p3 <- bias2plot(error38,
                title='model38 parameters bias',
                zrhos=c(expression(z(rho[12])),
                        expression(z(rho[34])),
                        expression(z(rho[14])),
                        expression(z(rho[23]))))
p4 <- bias2plot(error40,
                title='model40 parameters bias',
                zrhos=c(expression(z(rho[13])),
                        expression(z(rho[24])),
                        expression(z(rho[14])),
                        expression(z(rho[23]))))

(p1|p2)/(p3|p4)

CIF

beta1 <- c(-2,
           mean(coefs22$p1), mean(coefs36$p1),
           mean(coefs38$p1), mean(coefs40$p1))
beta2 <- c(-1.5,
           mean(coefs22$p2), mean(coefs36$p2),
           mean(coefs38$p2), mean(coefs40$p2))
gamma1 <- c(1.2,
            mean(coefs22$p3), mean(coefs36$p3),
            mean(coefs38$p3), mean(coefs40$p3))
gamma2 <- c(1,
            mean(coefs22$p4), mean(coefs36$p4),
            mean(coefs38$p4), mean(coefs40$p4))
w1 <- c(3,
        mean(coefs22$p5), mean(coefs36$p5),
        mean(coefs38$p5), mean(coefs40$p5))
w2 <- c(5,
        mean(coefs22$p6), mean(coefs36$p6),
        mean(coefs38$p6), mean(coefs40$p6))

t <- seq(30, 70.5, length.out=50)
delta <- 80
size <- length(beta1)
nt <- length(t)

cif1 <- vector(mode='list', length=size)
cif2 <- vector(mode='list', length=size)
for (i in seq(size)) {
    risklevel1 <- exp(beta1[i])
    risklevel2 <- exp(beta2[i])
    mlog1 <- risklevel1/(1+risklevel1+risklevel2)
    mlog2 <- risklevel2/(1+risklevel1+risklevel2)
    trajectory1 <- pnorm(w1[i]*atanh(2*t/delta-1)-gamma1[i])
    trajectory2 <- pnorm(w2[i]*atanh(2*t/delta-1)-gamma2[i])
    cif1[[i]] <- mlog1*trajectory1
    cif2[[i]] <- mlog2*trajectory2
}

model <- as_factor(c(rep('True', nt),
                     rep('model22', nt), rep('model36', nt),
                     rep('model38', nt), rep('model40', nt)))

data_cif <- tibble(curve=c(rep('Cause 1', 5*nt),
                           rep('Cause 2', 5*nt)), 
                   model=rep(model, 2), 
                   value=c(unlist(cif1), unlist(cif2)),
                   time=rep(t, 2*5))

ggplot(data_cif, aes(x=time, y=value, group=model))+
    geom_line(aes(linetype=model))+
    facet_wrap(~curve, labeller=label_bquote(CIF:.(curve)))+
    labs(x='Time', y=NULL, linetype='Model')+
    theme_minimal()+
    theme(axis.title.x=element_text(
              size=12,
              margin=unit(c(t=3, r=0, b=0, l=0), 'mm')),
          axis.text.x=element_text(size=11), 
          axis.text.y=element_text(size=11), 
          legend.position=c(0.11, 0.65),
          legend.text=element_text(size=12),
          legend.title=element_blank(), 
          legend.box.background=element_rect(color='black'), 
          strip.background=element_rect(colour='black', fill='white'),
          strip.text.x=element_text(size=12))

VCOV

plotS <- function(S, title) {
    longS <- reshape2::melt(S, na.rm=TRUE)
    longS$Var1 <- as_factor(longS$Var1)
    longS$Var2 <- fct_relevel(as_factor(longS$Var2), rev)
    ggplot(longS, aes(x=Var1, y=Var2, fill=value))+
        geom_tile(color='black', size=0.5)+
        geom_text(aes(label=round(value, 2)), size=4.5)+
        scale_fill_gradient(low='white', high='white')+
        labs(title=title)+
        scale_x_discrete(labels=c(expression(u[1]), 
                                  expression(u[2]), 
                                  expression(eta[1]),
                                  expression(eta[2])))+
        scale_y_discrete(labels=c(expression(eta[2]), 
                                  expression(eta[1]), 
                                  expression(u[2]),
                                  expression(u[1])))+
        theme_minimal()+
        theme(axis.text.x=element_text(size=12, color='black'),
              axis.text.y=element_text(size=12, color='black'), 
              axis.title.x=element_blank(), 
              axis.title.y=element_blank(),
              panel.grid.major=element_blank(),
              legend.position='none',
              plot.title=element_text(size=12, face='bold'))
}

p1 <- plotS(S=matrix(c(0.4, 0.15, 0.05, 0,
                       NA, 0.4, 0, 0.05,
                       NA, NA, 0.25, 0.1,
                       NA, NA, NA, 0.25), nrow=4, ncol=4),
            title='model22: true values')

p2 <- plotS(S=matrix(c(0.2, 0.15, 0.1, 0,
                       NA, 0.3, 0, 0.1,
                       NA, NA, 0.4, 0.15,
                       NA, NA, NA, 0.5), nrow=4, ncol=4),
            title='model22: initial guess')

means22 <- colMeans(coefs22)

p3 <- plotS(
    S=matrix(c(
        exp(means22[7]),
        tanh(means22[11])*sqrt(exp(means22[7]))*sqrt(exp(means22[8])),
        tanh(means22[13])*sqrt(exp(means22[7]))*sqrt(exp(means22[9])),
        0,
        NA, exp(means22[8]), 0,
        tanh(means22[14])*sqrt(exp(means22[8]))*sqrt(exp(means22[10])),
        NA, NA, exp(means22[9]),
        tanh(means22[12])*sqrt(exp(means22[9]))*sqrt(exp(means22[10])),
        NA, NA, NA, exp(means22[10])),
        nrow=4, ncol=4),
    title='model22: estimates')

p4 <- plotS(S=matrix(c(0.4, 0.15, 0.05, 0.2,
                       NA, 0.4, 0.2, 0.05,
                       NA, NA, 0.25, 0.1,
                       NA, NA, NA, 0.25), nrow=4, ncol=4),
            title='model36: true values')

p5 <- plotS(S=matrix(c(0.2, 0.15, 0.1, 0.2,
                       NA, 0.3, 0.2, 0.1,
                       NA, NA, 0.4, 0.15,
                       NA, NA, NA, 0.5), nrow=4, ncol=4),
            title='model36: initial guess')

means36 <- colMeans(coefs36)

p6 <- plotS(
    S=matrix(c(
        exp(means36[7]),
        tanh(means36[11])*sqrt(exp(means36[7]))*sqrt(exp(means36[8])),
        tanh(means36[13])*sqrt(exp(means36[7]))*sqrt(exp(means36[9])),
        tanh(means36[15])*sqrt(exp(means36[7]))*sqrt(exp(means36[10])),
        NA, exp(means36[8]),
        tanh(means36[16])*sqrt(exp(means36[8]))*sqrt(exp(means36[9])),
        tanh(means36[14])*sqrt(exp(means36[8]))*sqrt(exp(means36[10])),
        NA, NA, exp(means36[9]),
        tanh(means36[12])*sqrt(exp(means36[9]))*sqrt(exp(means36[10])),
        NA, NA, NA, exp(means36[10])), nrow=4, ncol=4),
    title='model36: estimates')

p7 <- plotS(S=matrix(c(0.4, 0.15, 0, 0.1,
                       NA, 0.4, 0.1, 0,
                       NA, NA, 0.25, 0.1,
                       NA, NA, NA, 0.25), nrow=4, ncol=4),
            title='model38: true values')

p8 <- plotS(S=matrix(c(0.2, 0.2, 0, 0.1,
                       NA, 0.3, 0.1, 0,
                       NA, NA, 0.4, 0.15,
                       NA, NA, NA, 0.5), nrow=4, ncol=4),
            title='model38: initial guess')

means38 <- colMeans(coefs38)

p9 <- plotS(
    S=matrix(c(
        exp(means38[7]),
        tanh(means38[11])*sqrt(exp(means38[7]))*sqrt(exp(means38[8])),
        0,
        tanh(means38[13])*sqrt(exp(means38[7]))*sqrt(exp(means38[10])),
        NA, exp(means38[8]),
        tanh(means38[14])*sqrt(exp(means38[8]))*sqrt(exp(means38[9])),
        0,
        NA, NA, exp(means38[9]),
        tanh(means38[12])*sqrt(exp(means38[9]))*sqrt(exp(means38[10])),
        NA, NA, NA, exp(means38[10])), nrow=4, ncol=4),
    title='model38: estimates')

p10 <- plotS(S=matrix(c(0.4, 0, 0.05, 0.2,
                        NA, 0.4, 0.2, 0.05,
                        NA, NA, 0.25, 0,
                        NA, NA, NA, 0.25), nrow=4, ncol=4),
             title='model40: true values')

p11 <- plotS(S=matrix(c(0.2, 0, 0.1, 0.2,
                        NA, 0.3, 0.2, 0.1,
                        NA, NA, 0.4, 0,
                        NA, NA, NA, 0.5), nrow=4, ncol=4),
             title='model40: initial guess')

means40 <- colMeans(coefs40)

p12 <- plotS(
    S=matrix(c(
        exp(means40[7]), 0,
        tanh(means40[11])*sqrt(exp(means40[7]))*sqrt(exp(means40[9])),
        tanh(means40[13])*sqrt(exp(means40[7]))*sqrt(exp(means40[10])),
        NA, exp(means40[8]),
        tanh(means40[14])*sqrt(exp(means40[8]))*sqrt(exp(means40[9])),
        tanh(means40[12])*sqrt(exp(means40[8]))*sqrt(exp(means40[10])),
        NA, NA, exp(means40[9]), 0,
        NA, NA, NA, exp(means40[10])), nrow=4, ncol=4),
    title='model40: estimates')

(p1|p2|p3)/(p4|p5|p6)/(p7|p8|p9)/(p10|p11|p12)
S <- matrix(c(1, 3, 5, 6,
              NA, 1, 6, 5,
              NA, NA, 2, 4,
              NA, NA, NA, 2), nrow=4, ncol=4)
longS <- reshape2::melt(S, na.rm=TRUE)
longS$Var1 <- as_factor(longS$Var1)
longS$Var2 <- fct_relevel(as_factor(longS$Var2), rev)
ggplot(longS, aes(x=Var1, y=Var2, fill=value))+
    geom_tile(color='black', size=0.5)+
    geom_text(aes( label=c('RISK\nLEVEL', 'RISK\nCORRELATION',
                           '1st ORDER\nRISK/TIME\nINTERACTION',
                           '2nd ORDER\nRISK/TIME\nINTERACTION',
                           'RISK\nLEVEL',
                           '2nd ORDER\nRISK/TIME\nINTERACTION',
                           '1st ORDER\nRISK/TIME\nINTERACTION',
                           'TRAJECTORY\nTIME',  'TIME\nCORRELATION',
                           'TRAJECTORY\nTIME')), size=5)+
    scale_fill_gradient(low='white', high='white')+
    scale_x_discrete(labels=c(expression(u[1]), 
                              expression(u[2]), 
                              expression(eta[1]),
                              expression(eta[2])))+
    scale_y_discrete(labels=c(expression(eta[2]), 
                              expression(eta[1]), 
                              expression(u[2]),
                              expression(u[1])))+
    theme_minimal()+
    theme(axis.text.x=element_text(size=15, color='black'),
          axis.text.y=element_text(size=15, color='black'), 
          axis.title.x=element_blank(), 
          axis.title.y=element_blank(),
          panel.grid.major=element_blank(),
          legend.position='none')